#ifndef SVM_KERNELS
#define SVM_KERNELS
/// Includes projects
#include "../common/framework.h"
#include "reduce.h"
/// Includes cuda
#include <vector_types.h>
#include <cuda_texture_types.h>
#include <texture_types.h>
#include <texture_fetch_functions.h>

texture<int2, cudaTextureType2D, cudaReadModeElementType> ellDataTex;
texture<int, cudaTextureType1D, cudaReadModeElementType> ellRowLenTex;

class Kernel {
public:
	__device__ float selfDot(int idxX);
	__device__ float dot(int idxX, float* Y);
	__device__ void parallelDot(unsigned int tid, float* pointerA, float* pointerAEnd, float* pointerB, float* sharedTemps);
	__device__ virtual ~Kernel() {}
};

class Linear : public Kernel {
public:
	__device__ float selfK(int idxX, float parameterA, float parameterB, float parameterC);
	__device__ float K(int idxX, float* Y, float parameterA, float parameterB, float parameterC, float* selfDot, int yIndex);
	__device__ void parallelKernel(unsigned int tid, float* pointerA, float* pointerAEnd, float* pointerB, float* sharedTemps, float parameterA, float parameterB, float parameterC);
	__device__ void dualK(int idxX, float *Y1, float *Y2, float parameterA, float parameterB, float parameterC, float &phiAB, float &phiAC, float* selfDot, int iHigh, int iLow);
};

class Polynomial : public Kernel {
public:
	__device__ float selfK(int idxX, float a, float r, float d);
	__device__ float K(int idxX, float* Y, float a, float r, float d, float* selfDot, int yIndex);
	__device__ void parallelKernel(unsigned int tid, float* pointerA, float* pointerAEnd, float* pointerB, float* sharedTemps, float a, float r, float d);
	__device__ void dualK(int idxX, float *Y1, float *Y2, float a, float r, float d, float &phiAB, float &phiAC, float* selfDot, int iHigh, int iLow);
};

class Gaussian : public Kernel {
public:
	__device__ float selfK(int idxX, float ngamma, float parameterB, float parameterC);
	__device__ float K(int idxX, float* Y, float ngamma, float paramterB, float parameterC, float* selfDot, int yIndex);
	__device__ void parallelKernel(unsigned int tid, float* pointerA, float* pointerAEnd, float* pointerB, float* sharedTemps, float ngamma, float parameterB, float parameterC);
	__device__ void dualK(int idxX, float *Y1, float *Y2, float ngamma, float parameterB, float parameterC, float &phiAB, float &phiAC, float* selfDot, int iHigh, int iLow);
};

class Sigmoid : public Kernel {
public:
	__device__ float selfK(int idxX, float a, float r, float parameterC);
	__device__ float K(int idxX, float* Y, float a, float r, float parameterC, float* selfDot, int yIndex);
	__device__ void parallelKernel(unsigned int tid, float* pointerA, float* pointerAEnd, float* pointerB, float* sharedTemps, float a, float r, float parameterC);
	__device__ void dualK(int idxX, float *Y1, float *Y2, float a, float r, float parameterC, float &phiAB, float &phiAC, float* selfDot, int iHigh, int iLow );
};

__device__ float Kernel::selfDot(int idxX) {
	float sum = 0.0f;
	for(int ptrX=0, lenX = tex1Dfetch(ellRowLenTex, idxX);ptrX<lenX;ptrX++) {
		int2 X = tex2D(ellDataTex, idxX, ptrX);
		float y = __int_as_float(X.y);
		sum += y*y;
	}
	return sum;
}

__device__ float Kernel::dot(int idxX, float* Y) {
	float sum = 0.0f;
	int lenX = tex1Dfetch(ellRowLenTex, idxX);
	for(int ptrX=0;ptrX<lenX;ptrX++) {
		int2 X = tex2D(ellDataTex, idxX, ptrX);
		sum += __int_as_float(X.y)* *(Y+X.x);
	}
	return sum;
}

__device__ void Kernel::parallelDot(unsigned int tid, float* pointerA, float* pointerAEnd, float* pointerB, float* sharedTemps) {
	pointerA += tid;
	pointerB += tid;
	sharedTemps[tid] = 0.0f;
	while (pointerA < pointerAEnd) {
		sharedTemps[tid] += (*pointerA) * (*pointerB);
		pointerA += blockDim.x;
		pointerB += blockDim.x;
	}
	__syncthreads();

	sumReduce<BLOCKSIZE>(sharedTemps,tid);
}

__device__ float Linear::selfK		(int idxX, float parameterA, float parameterB, float parameterC) {
	return selfDot(idxX);
}

__device__ float Polynomial::selfK	(int idxX, float a, float r, float d) {
	float sum = selfDot(idxX);
	sum = sum*a + r;
	return __powf(sum, d);
}

__device__ float Gaussian::selfK	(int idxX, float ngamma, float parameterB, float parameterC) {
	return 1.0f;
}

__device__ float Sigmoid::selfK		(int idxX, float a, float r, float parameterC) {
	float sum = selfDot(idxX);
	sum = sum*a + r;
	return tanhf(sum);
}

__device__ float Linear::K		(int idxX, float* Y, float parameterA, float parameterB, float parameterC, float* selfDot, int yIndex) {
	return dot(idxX,Y);
}

__device__ float Polynomial::K	(int idxX, float* Y, float a, float r, float d, float* selfDot, int yIndex) {
	float sum = dot(idxX,Y);
	sum = sum*a + r;
	return __powf(sum, d);
}

__device__ float Gaussian::K	(int idxX, float* Y, float ngamma, float paramterB, float parameterC, float* selfDot, int yIndex) {
	float sum = dot(idxX,Y);
	sum = -2.0f * sum;
	sum+= selfDot[idxX] + selfDot[yIndex];
	return __expf(sum*ngamma);
}

__device__ float Sigmoid::K		(int idxX, float* Y, float a, float r, float parameterC, float* selfDot, int yIndex) {
	float sum =  dot(idxX,Y);
	sum = sum*a + r;
	return tanhf(sum);
}

__device__ void Linear::parallelKernel		(unsigned int tid, float* pointerA, float* pointerAEnd, float* pointerB, float* sharedTemps, float parameterA, float parameterB, float parameterC) {
	parallelDot(tid, pointerA, pointerAEnd, pointerB, sharedTemps);
}

__device__ void Polynomial::parallelKernel	(unsigned int tid, float* pointerA, float* pointerAEnd, float* pointerB, float* sharedTemps, float a, float r, float d) {
	parallelDot(tid, pointerA, pointerAEnd, pointerB, sharedTemps);
	if (tid == 0) {
		float sum = sharedTemps[0]*a + r;
		sharedTemps[0] = __powf(sum, d);
	}
}

__device__ void Gaussian::parallelKernel	(unsigned int tid, float* pointerA, float* pointerAEnd, float* pointerB, float* sharedTemps, float ngamma, float parameterB, float parameterC) {
	pointerA += tid;
	pointerB += tid;
	sharedTemps[tid] = 0.0f;
	while (pointerA < pointerAEnd) {
		float diff = (*pointerA) - (*pointerB);
		sharedTemps[tid] += diff * diff;
		pointerA += blockDim.x;
		pointerB += blockDim.x;
	}
	__syncthreads();

	sumReduce<BLOCKSIZE>(sharedTemps,tid);
	if (tid == 0) {
		sharedTemps[0] = __expf(sharedTemps[0] * ngamma);
	}
}

__device__ void Sigmoid::parallelKernel		(unsigned int tid, float* pointerA, float* pointerAEnd, float* pointerB, float* sharedTemps, float a, float r, float parameterC) {
	parallelDot(tid, pointerA, pointerAEnd, pointerB, sharedTemps);
	if (tid == 0) {
		float sum = sharedTemps[0]*a + r;
		sharedTemps[0] = tanhf(sum);
	}
}

__device__ void Linear::dualK		(int idxX, float *Y1, float *Y2, float parameterA, float parameterB, float parameterC, float &phiAB, float &phiAC, float* selfDot, int iHigh, int iLow) {
	float sum1 = 0.0f;
	float sum2 = 0.0f;
	int lenX = tex1Dfetch(ellRowLenTex, idxX);
	for(int ptrX=0 ; ptrX<lenX ; ptrX++) {
		int2 X = tex2D(ellDataTex, idxX, ptrX);
		float fx = __int_as_float(X.y);
		sum1 += fx* *(Y1+X.x);
		sum2 += fx* *(Y2+X.x);
	}
	phiAB = sum1;
	phiAC = sum2;
}

__device__ void Polynomial::dualK	(int idxX, float *Y1, float *Y2, float a, float r, float d, float &phiAB, float &phiAC, float* selfDot, int iHigh, int iLow) {
	float sum1 = 0.0f;
	float sum2 = 0.0f;
	int lenX = tex1Dfetch(ellRowLenTex, idxX);
	for(int ptrX=0 ; ptrX<lenX ; ptrX++) {
		int2 X = tex2D(ellDataTex, idxX, ptrX);
		float fx = __int_as_float(X.y);
		sum1 += fx* *(Y1+X.x);
		sum2 += fx* *(Y2+X.x);
	}
	sum1 = sum1 * a + r;
	sum2 = sum2 * a + r;
	phiAB = __powf(sum1, d);
	phiAC = __powf(sum2, d);
}

__device__ void Gaussian::dualK		(int idxX, float *Y1, float *Y2, float ngamma, float parameterB, float parameterC, float &phiAB, float &phiAC, float* selfDot, int iHigh, int iLow) {
	float sum1 = 0.0f;
	float sum2 = 0.0f;

	int lenX = tex1Dfetch(ellRowLenTex, idxX);
	for(int ptrX=0 ; ptrX<lenX ; ptrX++) {
		int2 X = tex2D(ellDataTex, idxX, ptrX);
		float fx = __int_as_float(X.y);
		sum1 += fx* *(Y1+X.x);
		sum2 += fx* *(Y2+X.x);
	}
	float dotA = selfDot[idxX];
	{
		sum1 = -2.0f * sum1;
		float dotB = selfDot[iHigh];
		sum1 += dotA + dotB;
		phiAB = __expf(ngamma * sum1);
	}
	{
		sum2 = -2.0f * sum2;
		float dotC = selfDot[iLow];
		sum2 += dotA + dotC;
		phiAC = __expf(ngamma * sum2);
	}
}

__device__ void Sigmoid::dualK		(int idxX, float *Y1, float *Y2, float a, float r, float parameterC, float &phiAB, float &phiAC, float* selfDot, int iHigh, int iLow ) {
	float sum1 = 0.0f;
	float sum2 = 0.0f;
	int lenX = tex1Dfetch(ellRowLenTex, idxX);
	for(int ptrX=0 ; ptrX<lenX ; ptrX++) {
		int2 X = tex2D(ellDataTex, idxX, ptrX);
		float fx = __int_as_float(X.y);
		sum1 += fx* *(Y1+X.x);
		sum2 += fx* *(Y2+X.x);
	}
	sum1 = sum1 * a + r;
	sum2 = sum2 * a + r;
	phiAB = tanhf(sum1);
	phiAC = tanhf(sum2);
}

#endif
